Liquid-gas separation in colloidal electrolytes 
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The liquid-gas transition of an electroneutral mixture of oppositely charged colloids, studied 
by Monte Carlo simulations, is found in the low temperature - low density region. The critical 



temperature shows a non-monotonous behavior as a function of the interaction range, 



with a 



maximum at fca ~ 10, implying an island of coexistence in the K-p plane. The system is arranged 
in such a way that each particle is surrounded by shells of particles with alternating charge. In 
contrast with the electrolyte primitive model, both neutral and charged clusters are obtained in the 
vapor phase. 



PACS numbers: 82.70.Dd, 82.70.Gg, 64.70.Pf 

I. INTRODUCTION 

Electrostatic correlations play an important role in 
model, biological and applied systemgi, being their com- 
prehension a huge challenge for the liquid state re- 
searchers. While progress has been made in understand- 
ing the influence of those correlations in simpler ionic 
systems, only recently the importance of macromolec- 
ular correlations has been acknowledged. Furthermore, 
for charged colloidal systems, the interaction depends not 
only on the charge density and distribution, but also on 
the solvent properties; the interactions are thus tunable 
by acting on it (added salt, temperature...). This is the 
reason why phase diagrams are, in general, richer in col- 
loidal systems^. 

The Restricted Primitive Model (RPM) for ionic flu- 
ids is the simplest mixture showing strong effects due 
to charge correlations^. This model undergoes a liquid- 
gas transition at low temperature and density due 
to the strength of the correlationsii^ with Ising-like 
behaviors^,, despite the long range of the interactions. 
This model has been also extended to tackle size or charge 
asymmetriesSiLSiS, as an approach to charged colloids. 
The phase diagram of dipolar systems has been also stud- 
ied, where a vapor-liquid coexistence island is found as a 
function of the dipolar strengtbiS. 

In this work, we study the liquid-gas transition in the 
colloidal analogue of the RPM, using a mixture of op- 
positely charged colloids by means of computer simula- 
tions. We will focus on the case in which the concentra- 
tions of both colloidal species is the same, even though 
electroneutrality can be obeyed when the solution is not 
equimolar due to the electrolyte in the medium. Never- 
theless this is the case where the correlation effects be- 
tween unlike colloidal particles will be more significant. 
To date, the experimental works have focused on the 
crystal phases"'^^'"'^^, where superlattice crystals were re- 
ported, and recent simulations concentrated in the clus- 
tering of the particles^'^'^. In a previous work, the liquid- 



gas transition was indeed found for this system in the low 
density - low temperature regioni^, and we concentrate 
here in the effect of the range of the electrical interac- 
tion among particles, experimentally controlled by the 
electrolyte concentration in the medium. 

We find an unexpected non-monotonous behavior for 
the critical temperature with the range of interaction, 
with a maximum at kct w 10, thus showing a re-entrance 
phenomenon (fluid - phase separation - fluid) by increas- 
ing the salt concentration. This is a collective effect since 
a third term in the virial expansion of the pressure is 
necessary to describe the non-monotonous critical tem- 
perature curve. This fact is then rationalized from the 
strong correlation between oppositely charged particles, 
leading in both phases to particles surrounded by shells 
of alternating charges and clusters growing in the dilute 
phase. Under this configuration, repulsive and attractive 
interactions are screened in a different way, provoking a 
maximum in the critical temperature. Finally, we show 
that an increase of the amount of clusters is not a sig- 
nature of the proximity to the phase transition. That is 
actually driven by the energy gain in the dense phase to 
overcome the entropy lost in forming a liquid. 

The paper is organized as follows: Section II describes 
the details of the model and the simulation techniques. 
Section III is devoted to present and discuss the results 
of this work. Finally, the main conclusions can be found 
in Section IV. 



II. MODEL AND SIMULATIONS 

We simulate a 1:1 binary mixture of N spherical col- 
loidal particles, N /2 bearing a surface potential +(f) and 
N /2 with —(j). The interaction between colloids is mod- 
eled by the effective DLVO electrostatic interactiooi^, 
plus a hard-core repulsion. 



"^{^i]) = yHs{rij) + (l)t(f)je^p{-K{ri.j - cr)} (1) 



2 



where Vhs is the hard-sphere potential for a particle of 
diameter a, 0^ and (f>j are the surface potentials of the 
interacting particles (in appropriate units) and k is the 
inverse Debye length. This is obviously an approximation 
to the real problem, where the ions are not simulated, but 
an effective interaction between the colloids is used. This 
method, however, allows to study phase separation with 
large systems, as compared to those where the ions are 
explicitly includedi^. We shall use in this work reduced 
units: cr = 1, f7* = U/(j)'^, T* = keT/cl)'^ and the density 
p* — Na^ /V , where N is the number of particles, and V 
the volume of the system. 

Monte Carlo simulations have been used to compute 
the gas-liquid coexistence curve. First, Gibbs Ensemble 
Monte Carlo (GEMC) simulations with N = 432 parti- 
cles were run for different na values. Equilibration of the 
system takes very long times (1 — 5 • 10^ cycles); however, 
contrary to the RPM, multiparticle moves do not speed 
up the the equilibration rate due to the high density of 
the liquid phase. Production runs comprised at least 10^ 
cycles, taken at equilibrium. 

Grand Canonical Monte Carlo (GCMC) simulations 
aided with reweighting techniquesi^iii were used to ver- 
ify the results from GEMC. Owing to the long range of 
the correlations, the coexistence curve depends strongly 
on the system size in the neighborhood of the critical 
point. This fact can be used to locate the critical point 
using the mixed-field finite size scaling analysis devel- 
oped by Bruce and Wilding^-. This method is based on 
the asymmetry of the density distribution, reflecting the 
absence of particle-hole symmetry in off-lattice models. 
Thus, the order parameter of the transition is not the 
density, but a mixture between the density and the en- 
ergy: M ~ yO -|- SM (where s is a system-dependent field 
mixing parameter). Precisely at criticality, the distribu- 
tion of M in a large enough system, with box size L, 
takes on the universal form: 

Pl{M) ^ aLP* {aL[M - {M)]) (2) 

where P* (x) is an universal distribution function for each 
universality class, < ... > is the average evaluated at 
critical conditions and a/, ~ L^^'^ . Thus, using GCMC 
simulations, critical parameters for different box sizes are 
calculated, and applying the corresponding scaling laws, 
the critical parameters for an infinite system can be ob- 
tained 

T,(L) - r,(oo) ^ L-^^+i)/'^ (3) 
p,{L) ~ p,(oo) ^ (4) 

where u and a are the critical exponents associated to the 
correlation length and heat capacity divergence, respec- 
tively, and 9 is the correction-to-scaling exponent. For 
the 3D Ising universality class, a ~ 0.11, v « 0.629 and 
« 0.54i^. Some remarks are pertinent at this point. 
First of all, the Bruce- Wilding method does not identify 
the universality class of the system but to obtain the criti- 
cal parameters assuming certain universality class. Since 




0.131 — ' — I — ' — I — ' — I — ' — I — ' — I — ' — I — ' — I 
0.1 0.2 0.3 0.4 0.5 0.6 

P* 

FIG. 1: Gas-liquid coexistence points for different values of 
K = 3.9 circles, k, = Q, open squares, k = 10, upward triangles, 
K = 15, downward triangles, and k = 20, diamonds. The 
points without error bars mark the critical points, estimated 
with the law of rectilinear diameters. 



our model presents a short range interaction, 3D Ising 
criticality is expected. On the other hand, the Bruce- 
Wilding method is not consistent with the Yang- Yang 
anomaly. This method can be improved including the 
pressure as a scaling field, but for Ising-type systems this 
contribution is usually negligible^". Therefore, we follow 
the Bruce- Wilding method. 

To carry out this analysis, we have used GCMC sim- 
ulations with L = 8, 9, 10 and 12 system sizes. Runs 
comprised 25 • 10^ configurations for L = 8 and L = 9 
and 75 • 10^ configurations for L — IQ and L — 12, con- 
figurations are spaced by 2 < > attempts to insert or 
remove a random particle. Statistical errors were taken 
from three independent simulations runs. 



III. RESULTS AND DISCUSSION 

The liquid-gas transition, studied by means of Gibbs 
Ensemble Monte Carlo and Grand Canonical simulations, 
is found in the low-T-low-p region, where strong correla- 
tions between oppositely charged particles emerge, pro- 
ducing a thermodynamic instabilitjsiii^. We find liquids 
and vapors composed of the same number of positive 
and negative particles on average. The GEMC results 
for different k are presented in Fig. In contrast to 
monocomponent attractive systems^, the critical tem- 
perature evolves non-monotonically. At low k, demixing 
occurs at higher temperatures the shorter the interac- 
tion range, and only at high k, the critical temperature 
moves to lower values as k increases. Noteworthily, this 
behavior implies a closed region of phase separation in 
the constant-T, n-p plane, which is experimentally more 
accessible than the T-p one for colloids. A reentrance 
phenomenon is thus predicted by increasing the salt con- 
centration, from fluid to phase separation and to fluid 
again, similarly to the dipolar systemsiSi. 
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To verify the dependence of the critical temperature 
with the salt concentration, we have also used GCMC 
simulations what are more precise than the GEMC ones. 
Fig. 121 (upper panel) shows Pl (M) matched to the uni- 
versal 3D Ising distribution for kct = 6 with different 
box sizes: L = 8,9, 10 and 12. This matching gives 
us a critical temperature for each system size, that can 
be used to obtain Tc in an infinite system as the lin- 
eal fit to the suitable scaling law^^ (results presented 
in the lower panel of Fig. I^J. Note the small de- 
pendence of the critical temperature with the size of 
the simulation boxes; the critical density also depends 
slightly on L (not shown). These simulations confirm 
the behavior of the critical temperature with the salt 
concentration: Tc increases when k increases at low k 
{Tc{Ka = 6) < Tc{k<t — 10)) and decreases at high 
enough salt concentration {Tc{k<t — 15) < Tc{k(t — 10)). 

The size dependence of the non universal parameter 
ul L^^'^ (not shown) permits an estimation of the 
ratio P/iy. The values obtained are similar to those of the 
Ising model (fi/f = 0.518) for the three ranges studied: 
P/v = 0.523(14), 0.509(23) and 0.503(20) for na = 6,10 
and 15, respectively. The collapse of the data onto the 
3D-Ising curve (solid line) in Fig. |21 the agreement of 
numerical values for the ratios /3/iy, and the analysis of 
ttL, support the compatibility with 3D-Ising criticality, 
as expected for short-range potentials. 

The critical temperature and density from both GEMC 
and GCMC simulations are presented in Fig. |21 as a 
function of k. The agreement between both simula- 
tion results confirms the behavior described above; on 
the other hand, the critical density, which is better es- 
timated with the GCMC simulations, increases with k. 
The critical parameters have been calculated using the 
virial expansion of the pressure up to second and third 
orders (continuous and dashed lines in Fig. 13 respec- 
tively): /?P(p) = (3Pgl{p) + B2p\+ B^p^ where PgUp) 
is the Carnahan-Starling expression for the pressure of 
hard spheres^^i^. The second order expansion produces 
a liquid-gas transition, driven by the attractions between 
oppositely charged particles, which behaves monotoni- 
cally with kct (solid lines). On the other hand, the ex- 
pansion up to third order correctly reproduces the maxi- 
mum in Tc vs. Ka (dashed lines). Because is needed to 
reproduce qualitatively the results, the increasing trend 
of Tc vs. Ka at long interaction ranges comes from the 
interaction between, at least, three particles, mainly the 
repulsion between two particles with similar charge at- 
tracted to a third one. This repulsion is more important 
for lower k, whereas the attraction energy is not greatly 
increased. The origin of the maximum is, thus, similar 
to that of the dipolar systemiS. These results also agree 
qualitatively with previous calculations for an attractive- 
repulsive mixture of Yukawa potentials using an expan- 
sion of the internal energy and the equation of state, 
based on the mean-spherical approximation (MSA)2&. 

As in the ionic fluids, the demixing in a dense and a 
dilute phase is driven by the strong correlations between 
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FIG. 2: Upper panel: order parameter probability distribu- 
tion Pl{x) for Kcr = 6 and several box sizes: L = 8 (black 
circles), L = 9 (red squares) L = 10 (green diamonds) and 
L = 12 (blue triangles); solid line for 3D-Ising model. Lower 
panel: variation of the critical temperature with the box sizes 
for KCT = 6 (black squares) kct = 10 (red squares) and no — 15 
(green diamonds). Lines are lineal fits to the numerical re- 
sults. 
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FIG. 3: Critical temperature and density (inset) as a function 
of Ka: solid squares from GEMC and open red diamonds from 
GCMC simulations. Black lines are the results using the virial 
expansion for the pressure up to second (solid line) and third 
order (broken line) in density. Green lines are results from 
MSA on a Yukawa mixture--. 
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FIG. 4: Partial pair distribution functions for T = 0.18 and 
p = 0.6 with K = 2 (black curve), k = Q (red curve) and 
ft = 15 (green curve); opposite sign (main figure) and same 
sign (inset). 



oppositely charged particle&WA, \Ye have investigated 
the internal correlations in two supercritical {T* — 0.18) 
states, p* = 0.6 and p* = 0.05, marked in Fig. ^as as- 
terisks. Fig. 0] plots the opposite sign, (r), and same 

sign (inset), g++{r) = g — (r), contributions to the pair 
distribution function for the dense state (with kct = 2, 6 
and 15). The system is composed of layers with alternat- 
ing sign particles surrounding every particle; first layer of 
particles of opposite sign, second one of the same sign... 
This structure extends to distances of several radii, much 
larger than the interaction range. As observed in the Fig- 
ure, the peaks shift to shorter distances as k increases. 
Qualitatively similar features have been predicted using 
the MSA for the partial structure factors^. 

Second neighbors of one particle comprise the first 
shell of particles with the same sign, which moves to 
larger distances as k decreases. This is due to the repul- 
sion between similarly charged particles and correlations 
through unlike particles. The first neighbor layer, thus, 
drive the gas liquid transition, whereas the second one 
impedes it, resulting in a maximum of the critical tem- 
perature. A peak at r = 2(t marks the presence of linear 
arrangements (three particles long) at all ranges studied. 
This peak is not caused by the repulsion of the central 
particle (because the repulsion range is too short), but 
due to the repulsion of the particles in the second layer, 
first peak in g++{r). Visual inspection of the system 
shows that these strings are distributed randomly (thus 
cannot attributed to crystallites), and are only three par- 
ticles long (no peak a,t r — 3a in g{r)). 

Now we move to the dilute supercritical state T* — 
0.18 and p* = 0.05, presented in Fig. Ofor different in- 
teraction ranges (upper panel). The system is composed 
of clusters of particleai^, and the fraction of particles in 
every cluster, x„, is presented in the lower panel of the 
figure for different values of k. Again the layering of 
particles inside clusters is observed although only two or 
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FIG. 5: Upper panel: Partial pair distribution functions for 
T* = 0.18 and p = 0.05 with k. ^ 2 (black curve), k = 6 
(red curve), k = 15 (green curve) and it — 30 (blue curve); 
oppositely sign (main figure) and same sign (inset). Lower 
panel: Fraction of particles in a cluster with n particles for 
the same state and different k, as labeled. Xn is normalized 
as: J2^ri = 1- 



three layers are seen due to the finite size of the clus- 
ters. The layers move closer as k increases, the effect 
being more dramatic than in the dense state due to the 
lower density of the system. Notably, the absence of the 
peak at r = 2cr confirms the collective origin of the linear 
structures found in the dense phase (Fig. 

Different from the RPM22i2S, the distribution of clus- 
ters is monotonous, containing both neutral and charged 
clusters, no matter how close to the transition. However, 
this distribution shows also a non-monotonic behavior 
with k; for long ranges the number of large clusters in- 
creases when K increases, whereas the opposite trend is 
observed at high k, the maximum number of big clusters 
found for kct = 6.0 (for not-too-big clusters, the maxi- 
mum is at Ka = 3.9). These values do not agree with the 
maximum of T*, which indicates that maximal proxim- 
ity to the transition does not imply maximal tendency to 
form large clusters in dilute systems. 

Finally, to complete the understanding of the phase 
transition and its driving mechanism. Fig. |S1 depicts the 
internal (electrostatic) energy with its contributions from 
attractions and repulsions. The same supercritical states 
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FIG. 6: Electrostatic energy (squares), and contributions 
from repulsions (upwards triangles) and attractions (down- 
wards triangles) for two states as labeled. 

are presented (T* = 0.18; p* = 0.6 (upper panel) and 
p* = 0.05 (lower panel)). At low k, the total energy com- 
prises repulsive and attractive contributions. For larger 
K, however, the energy contains only the attractive con- 
tributions due to the different distances between similar 
sign and opposite sign pairs (see Fig. Therefore, the 
behavior at high k is similar to that of a monocompo- 
nent system. For the dilute state, the clustering implies 
a low number of bonds, as compared to the dense state, 
resulting in a lower energy (in absolute value) . The mini- 
mum in the energy can be rationalized using simple small 
clusters, i.e. trimers and tetramers^S. 

Comparison of the energy curves with the behavior of 
critical temperature shows that the driving mechanism 



for liquid-gas separation in this system is the energy gain 
in forming a dense (liquid) phase, and not the ability to 
form large clusters in dilute phases (see Fig. EI), as pro- 
posed for the asymmetric ionic fluids'^^. The effect arises 
from the difference in the number of interacting pairs of 
similarly charged colloids in the dilute and dense phases. 
The internal energy calculated using the MSA for a mix- 
ture of Yukawa potentials shows similar trends, with min- 
ima moving to higher k as the density increases^^. 



IV. CONCLUSIONS 

In conclusion, we have studied the liquid-gas transi- 
tion in a system which can be considered as the colloidal 
analogue of the ionic fluid, i.e. a 1:1 mixture of oppo- 
sitely charged colloids. The coexistence region is found 
in the low density-low temperature region, with 3D-Ising 
criticality. The critical temperature shows a non mono- 
tonic behavior with the range of the interaction, k; in- 
creasing in low values of k and decreasing for the higher 
ones. This prediction implies a closed region of liquid-gas 
demixing at constant temperature in the K — p plane. The 
condensed phases are arranged in such a way that each 
particle is surrounded by shells of particles with alternat- 
ing charge. The internal energy shows that the overall 
trend of the critical temperature is lead by the energy 
gain in forming a dense phase. Liquid-gas separation 
at long ranges is hindered by the repulsion between simi- 
larly charged particles bonded to a third one (of opposite 
sign), whereas at high k this repulsion is negligible. 
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